package app.jppfnetmorph;

import static optimization.LinearAlgebra.*;

import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStreamReader;
import java.util.ArrayList;
import java.util.List;
import java.util.StringTokenizer;

import neuron.NetworkStatisticsProxy;
import optimization.Logger;
import theory.vanpelt.bifurcation.BEFindValues;
import optimization.LinearAlgebra;
import optimization.RetrospectiveApproximation;
import netmorph.model.parts.ObliqueElongPart;
import netmorph.model.parts.TuftBEPart;
import netmorph.model.parts.BasalBEPart;
import netmorph.model.parts.BasalElongPart;
import netmorph.model.parts.NetmorphModelBESPart;
import netmorph.model.parts.NetmorphModelElongPart;
import netmorph.model.parts.ObliqueBEPart;
import netmorph.model.parts.TuftElongPart;

import org.jppf.client.JPPFJob;

import netmorph.model.NetmorphModel;
import netmorph.model.parts.BESInputModel;
import netmorph.model.parts.ElongInputModel;

import common.notification.EmailNotifier;
import common.statistics.Series;

public class RAConsole {

	double[] h = new double[]{ 0.1, 0.05};

	RetrospectiveApproximation ra;
	NelderMeadParallel nmp;
	
	String name;

	int MAX_ITER_RA = 10;
	
	private double[] yt;

	private char mode = 'b';
	private NetworkStatisticsProxy.Selector sel = NetworkStatisticsProxy.basalTS;
	private Class<? extends NetmorphModelBESPart> partClass = BasalBEPart.class;
	private Class<? extends NetmorphModelElongPart> elongClass = BasalElongPart.class;
	
	private double[] basalBE = new double[] { 1.72, 0.42 };
	private double[] tuftBE = new double[] { 6.0, 0.25 };
	private double[] obliqueBE = new double[] { 1.24, 0.32 };
	
	private double[] basalElong = new double[]   { 0.000042593, 0.000042593, 10, 10 };
	private double[] tuftElong = new double[]  {   0.000346296, 0.000037037, 10, 10 };
	private double[] obliqueElong = new double[] { 0.000051852, 0.000046296, 10, 10 };
	
	private NetmorphConfig.Initializer init = new NetmorphConfig.Initializer() {
		@Override
		public void init(NetmorphModel model, BESInputModel bes, ElongInputModel elong) {
			
			model.setParam("pyramidal.min_basal", "7");
			model.setParam("pyramidal.max_basal", "7");

			model.setParam("APD.tsem.obliques.PDF.mean", "7");
			model.setParam("APD.tsem.obliques.PDF.std", "0");

			
			bes.setBasalB(basalBE[0]);
			bes.setBasalE(basalBE[1]);

			bes.setTuftB(tuftBE[0]);
			bes.setTuftE(tuftBE[1]);
			
			bes.setObliqueB(obliqueBE[0]);
			bes.setObliqueE(obliqueBE[1]);
			
			elong.setBERIPDFMean(basalElong[0]);
			elong.setBERIPDFStdev(basalElong[1]);
			elong.setBTSEMBPDFMean(basalElong[2]);
			elong.setBTSEMBPDFStdev(basalElong[3]);

			elong.setTERIPDFMean(tuftElong[0]);
			elong.setTERIPDFStdev(tuftElong[1]);
			elong.setTTSEMBPDFMean(tuftElong[2]);
			elong.setTTSEMBPDFStdev(tuftElong[3]);

			elong.setOERIPDFMean(obliqueElong[0]);
			elong.setOERIPDFStdev(obliqueElong[1]);
			elong.setOTSEMBPDFMean(obliqueElong[2]);
			elong.setOTSEMBPDFStdev(obliqueElong[3]);
		}
	};
	
	VectorFunction beCons = new VectorFunction() {
		@Override
		public int yDim() {
			return 2;
		}
		@Override
		public double[] evaluate(double[] p) {
			if (p[0] < 1) p[0] = 1;
			if (p[0] > 6) p[0] = 6;
			if (p[1] < 0) p[1] = 0;
			if (p[1] > 1) p[1] = 1;
			return p;
		}
	};

	private void setTarget(double[] yt)
	{
		this.yt = yt;
		
		// init RA
		double[] xinit = BEFindValues.find(yt[0], yt[1]);
		ra = new RetrospectiveApproximation(null, beCons, xinit, h) {
			protected int sampleSize(int iter)
			{
				return (int)(200*Math.pow(2,iter)); 
			}
		};
		
	}

	private void setInit(double[] x)
	{
		switch (mode) {
			case 'b': basalBE = x; break;
			case 't': tuftBE = x; break;
			case 'o': obliqueBE = x; break;
		}
		System.out.println("INIT " +
				"basal " + vPrint(basalBE) +
				"tuft " + vPrint(tuftBE) +
				"oblique " + vPrint(obliqueBE) 
				);
	}
	
	private void step() throws Exception
	{
		ra.beginIteration();
		nmp = new NelderMeadParallel(sel, partClass, init, ra.N, yt, ra.term, beCons);
		ra.xa = nmp.search(ra.xa, vMul(ra.hscale, h), 10);
		Logger.format("[%s] RA%d   X %.5f, %.5f   err %.5f   (N=%d, rnd=%d, h=%.5f, nu=%.5f)", name, ra.iter, ra.xa[0], ra.xa[1], nmp.getEstimateError(), ra.N, ra.randSeed, ra.hscale, ra.term);
		ra.endIteration();
		EmailNotifier.notifyNoFail(Logger.getLog());
	}

	private String eval(double x, double y, int n) throws Exception 
	{
		List<double[]> points = new ArrayList<double[]>();
		double[] p = new double[]{ x, y }; 
		points.add(p);
		JPPFJob job = NelderMeadParallel.compJob(points, n, "ConsoleEval", partClass, init);
		List<JobResult> jrl = JPPFRunner.runJob(job);
		Series res = NelderMeadParallel.collect(jrl, p, sel);

		EmailNotifier.notifyNoFail(String.format("f(%.5f, %.5f, n=%d) -> %s\n", x, y, n, res.formatNormDist()));

		Series ts = null, is = null;
		switch(mode) {
			case 'b':
				ts = NelderMeadParallel.collect(jrl, p, NetworkStatisticsProxy.basalTSLen);
				is = NelderMeadParallel.collect(jrl, p, NetworkStatisticsProxy.basalISLen);
				break;
			case 't':
				ts = NelderMeadParallel.collect(jrl, p, NetworkStatisticsProxy.tuftTSLen);
				is = NelderMeadParallel.collect(jrl, p, NetworkStatisticsProxy.tuftISLen);
				break;
			case 'o':
				ts = NelderMeadParallel.collect(jrl, p, NetworkStatisticsProxy.obliqueTSLen);
				is = NelderMeadParallel.collect(jrl, p, NetworkStatisticsProxy.obliqueISLen);
				break;
		}			
		StringBuilder sb = new StringBuilder();
		sb.append("termseg = [\n");
		for (double d : ts.values()) { sb.append(String.format("%.4f ", d)); }
		sb.append("];\n\n");
		sb.append("intseg = [\n");
		for (double d : is.values()) { sb.append(String.format("%.4f ", d)); }
		sb.append("];\n\n");
		EmailNotifier.notifyNoFail(sb.toString());

		return res.formatNormDist();
	}

	private double[] compElong(double[] p, int n) throws Exception 
	{
		List<double[]> points = new ArrayList<double[]>();
		points.add(p);
		
		Class<? extends NetmorphModelElongPart> part = null;
		if (mode == 'b') part = BasalElongPart.class;
		if (mode == 't') part = TuftElongPart.class;
		if (mode == 'o') part = ObliqueElongPart.class;
		JPPFJob job = NelderMeadParallel.compJob(points, n, "Elong", null, part, init);			
		List<JobResult> jrl = JPPFRunner.runJob(job);

		Series ts = null, is = null;
		switch(mode) {
			case 'b':
				ts = NelderMeadParallel.collect(jrl, p, NetworkStatisticsProxy.basalTSLen);
				is = NelderMeadParallel.collect(jrl, p, NetworkStatisticsProxy.basalISLen);
				break;
			case 't':
				ts = NelderMeadParallel.collect(jrl, p, NetworkStatisticsProxy.tuftTSLen);
				is = NelderMeadParallel.collect(jrl, p, NetworkStatisticsProxy.tuftISLen);
				break;
			case 'o':
				ts = NelderMeadParallel.collect(jrl, p, NetworkStatisticsProxy.obliqueTSLen);
				is = NelderMeadParallel.collect(jrl, p, NetworkStatisticsProxy.obliqueISLen);
				break;
		}			
		double[] y = new double[]{ ts.mean(), ts.stdev(), is.mean(), is.stdev() };
		return y;
	}
	
	
	private void elongjacob(int n) throws Exception 
	{
		double[] p = null;
		if (mode == 'b') p = basalElong.clone(); 
		if (mode == 't') p = tuftElong.clone(); 
		if (mode == 'o') p = obliqueElong.clone(); 
		if (p == null) return;
		
		double[] y0 = compElong(p, n);
		double[] w = { 10000, 10000, 1, 1 };
		
		for (int k = 0; k < 4; k++) {
			double[] p1 = p.clone();
			p1[k] *= 2;
			double[] y1 = compElong(p1, n);

			double[] p2 = p.clone();
			p2[k] /= 2;
			double[] y2 = compElong(p2, n);
			
			double[] dy = vSub(y1, y2);
			double dp = (p1[k] - p2[k])*w[k];

			System.out.println("dy : " + LinearAlgebra.vPrintL(dy) + ", dp : " + dp);
			System.out.println("dy/dp: " + LinearAlgebra.vPrintL(vMul(1/dp, dy)));
		}
		
	}
	
	private void searchElong3(double[] yt) throws Exception {

		double[] p = null;
		if (mode == 'b') p = basalElong.clone(); 
		if (mode == 't') p = tuftElong.clone(); 
		if (mode == 'o') p = obliqueElong.clone(); 

		int n = 200;
		String xx = null;
		
		for (int i = 0; i < 5; i++) {
			
			double[] w = { 0.5, 0.5, 0.5, 0.5 };

			for (int k = 0; k < 4; k++) { 

				for (int j = 0; j < 10; j++) {

					System.out.println(k+"."+j);
					
					double[] y0 = compElong(p, n);
					double err0 = vNorm2(vSub(y0, yt));

					double[] x1 = p.clone();
					x1[k] += w[k]*p[k]; 
					double[] y1 = compElong(x1, n);
					double err1 = vNorm2(vSub(y1, yt));
	
					double[] x2 = p.clone();
					x2[k] -= w[k]*p[k]; 
					double[] y2 = compElong(x2, n);
					double err2 = vNorm2(vSub(y2, yt));
	
					if (err1 < err2 && err1 < err0) { p = x1; } 
					else if (err2 < err1 && err2 < err0) { p = x2; } 
					else { w[k] /= 2; }
				}
			}
//			w = vMul(0.5, w);

			double[] yk = compElong(p, n);
			
			xx = String.format("[%s] 2EL%2d  ", name, i);
			xx += " X " + vPrintL(p) + ",  Y  " + vPrint(yk) + ",  TARGET " + vPrint(yt);
			System.out.println(xx);
			

//			n = n * 2;
		}
		EmailNotifier.notifyNoFail(xx);
	}

	/**
	 * 
	 * @param yt1  tsm, tss, ism, iss
	 * @throws Exception
	 */
	private void search2(double[] yt1) throws Exception {

		double[] p = null;
		if (mode == 'b') p = basalElong.clone(); 
		if (mode == 't') p = tuftElong.clone(); 
		if (mode == 'o') p = obliqueElong.clone(); 
		double[] eh = { 0.1, 0.1, 0.1, 0.1 };
		double[] t = { 1, 1, 1, 1 };
		int n = 200;
		
		for (int i = 0; i < 20; i++) {

			double[] yb = compElong(p, n);
			double errb = vNorm2(vSub(yb, yt1));
			
			for (int k = 0; k < 4; k++) { 

				double[] x1 = vMul((1+eh[k]), p); 
				double[] y1 = compElong(x1, n);
				double err1 = vNorm2(vSub(y1, yt1));

				double[] x2 = vMul((1-eh[k]), p); 
				double[] y2 = compElong(x2, n);
				double err2 = vNorm2(vSub(y2, yt1));

				if (err1 < err2 && err1 < errb) {
					// accept 1
					p = x1;
				}
				else if (err2 < err1 && err2 < errb) {
					// accept 2
					p = x2;
				}
				else {
					// reject, halve h
					eh[k] = eh[k] / 2;
				}
			}			

			String xx = String.format("[%s] EL%2d  ", name, i);
			xx += " X " + vPrintL(p) + ",  Y  " + vPrint(yb) + ",  TARGET " + vPrint(yt1);
			xx += ";; H " + vPrintL(eh) + ", T " + vPrintL(t);
			System.out.println(xx);
			
			if (i % 10 == 9)
				EmailNotifier.notifyNoFail(xx);

//			n = n * 2;
		}		
	}
	
	private void search() throws Exception {

		while(true) {

			System.out.println("SEARCH STEP");
			step();

			System.out.println("SEARCH ADV");
			if (ra.iter >= MAX_ITER_RA) {
				double[] xest = ra.estimateRoot();
				Logger.format("[%s] RA--   X %.5f, %.5f   ", name, xest[0], xest[1]);
				EmailNotifier.notifyNoFail(Logger.getLog());
				return;
			}
		}
	}
	

	public void run() throws IOException
	{
		System.out.println("RA Console v1.00 - 09-08-2010");
		System.out.println("Commands: target x y, search, name xx, eval x y n, mode basal|oblique|tuft, init x y, rainit x y, exit");
		String line = "";
		BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
		while ( (line = br.readLine()) != null) {
			StringTokenizer st = new StringTokenizer(line.trim().toLowerCase());
			String cmd = st.nextToken();
			if (cmd.equals("exit")) {
				System.exit(1);
			}
			else if (cmd.equals("target")) {
				double x = Double.parseDouble(st.nextToken());
				double y = Double.parseDouble(st.nextToken());
				System.out.println("TARGET " + x + ", " + y);
				setTarget(new double[] { x, y });
			}
			else if (cmd.equals("rainit")) {
				double x = Double.parseDouble(st.nextToken());
				double y = Double.parseDouble(st.nextToken());
				System.out.println("RAINIT " + x + ", " + y);
				ra.xa = new double[] { x, y };
			}
			else if (cmd.equals("name")) {
				name = st.nextToken();
				System.out.println("SET NAME = " + name);
			}
			else if (cmd.equals("elongjacob")) {
				int n = Integer.parseInt(st.nextToken());
				try { elongjacob(n); } 
				catch (Exception e) {
					e.printStackTrace();
				}
			}
			else if (cmd.equals("mode")) {
				String tmp = st.nextToken();
				if (tmp.equalsIgnoreCase("basal")) {
					sel = NetworkStatisticsProxy.basalTS;
					partClass = BasalBEPart.class;
					mode = 'b';
				}
				else if (tmp.equalsIgnoreCase("tuft")) {
					sel = NetworkStatisticsProxy.tuftTS;
					partClass = TuftBEPart.class;
					mode = 't';
				}
				else if (tmp.equalsIgnoreCase("oblique")) {
					sel = NetworkStatisticsProxy.obliqueTS;
					partClass = ObliqueBEPart.class;
					mode = 'o';
				}
				System.out.println("SELECTOR = " + sel);
				System.out.println("MODELPART = " + partClass);
			}
			else if (cmd.equals("init")) {
				double x = Double.parseDouble(st.nextToken());
				double y = Double.parseDouble(st.nextToken());
				System.out.println("INIT " + x + ", " + y);
				setInit(new double[] { x, y });
			}
			else if (cmd.equals("initall")) {
				double a = Double.parseDouble(st.nextToken());
				double b = Double.parseDouble(st.nextToken());
				double c = Double.parseDouble(st.nextToken());
				double d = Double.parseDouble(st.nextToken());
				double e = Double.parseDouble(st.nextToken());
				double f = Double.parseDouble(st.nextToken());
				System.out.println("INIT-ALL " + a +"," +b+","+c+","+d+","+e+","+f);
				basalBE = new double[]{a, b};
				tuftBE = new double[]{c, d};
				obliqueBE = new double[]{e, f};
			}
			else if (cmd.equals("search")) {
				System.out.println("SEARCH");
				try {
					search();
				} catch (Exception e) {
					e.printStackTrace();
				}
			}
			else if (cmd.equals("searchelong")) {
				double x = Double.parseDouble(st.nextToken());
				double y = Double.parseDouble(st.nextToken());
				double z = Double.parseDouble(st.nextToken());
				double w = Double.parseDouble(st.nextToken());
				System.out.println("SEARCH-ELONG " + x+ ", " + y +", " + z + ", " + w);
				try {
					search2(new double[]{ x, y, z, w });
				} catch (Exception e) {
					e.printStackTrace();
				}
			}
			else if (cmd.equals("s3")) {
				double x = Double.parseDouble(st.nextToken());
				double y = Double.parseDouble(st.nextToken());
				double z = Double.parseDouble(st.nextToken());
				double w = Double.parseDouble(st.nextToken());
				System.out.println("SEARCH-ELONG-2 " + x+ ", " + y +", " + z + ", " + w);
				try {
					searchElong3(new double[]{ x, y, z, w });
				} catch (Exception e) {
					e.printStackTrace();
				}
			}
			else if (cmd.equals("eval")) {
				double x = Double.parseDouble(st.nextToken());
				double y = Double.parseDouble(st.nextToken());
				int n = Integer.parseInt(st.nextToken());
				System.out.println("EVAL " + x + ", " + y + " (N=" + n + ")");
				try {
					String res = eval(x,y,n);
					System.out.println("F(" + x +"," + y + ") = " + res+ " (N=" + n + ")");
				} catch (Exception e) {
					e.printStackTrace();
				}
			}
			else {
				System.out.println("SYNTAX ERROR: " + line.trim());
			}
		}
	}
	
	
	public RAConsole() {
	}


	
	public static void main(String[] args) throws IOException {
		System.setProperty("jppf.config", "config/jppf.properties");
		System.setProperty("log4j.configuration", "config/log4j.properties");
		RAConsole rac = new RAConsole();
		rac.run();
	}
}
